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Abstract 

Large-scale clustering of highly biased tracers of large-scale structure has emerged as one of the 
best observational probes of primordial non-Gaussianity of the local type (i.e. /j^f 1 ). This type 
of non-Gaussianity can be generated in multifield models of inflation such as the curvaton model. 
Recently, Tseliakhovich, Hirata, and Slosar showed that the clustering statistics depend qualitatively 
on the ratio of inflaton to curvaton power £ after reheating, a free parameter of the model. If £ is 
significantly different from zero, so that the inflaton makes a non- negligible contribution to the 
primordial adiabatic curvature, then the peak-background split ansatz predicts that the halo bias 
will be stochastic on large scales. In this paper, we test this prediction in iV-body simulations. We 
find that large-scale stochasticity is generated, in qualitative agreement with the prediction, but 
that the level of stochasticity is overpredicted by ~30%. Other predictions, such as £ independence 
of the halo bias, are confirmed by the simulations. Surprisingly, even in the Gaussian case we do 
not find that halo model predictions for stochasticity agree consistently with simulations, suggesting 
that semi-analytic modeling of stochasticity is generally more difficult than modeling halo bias. 



1 Introduction 



One of the most exciting prospects for cosmology in the near future is the ability to constrain 
the physics of inflation [1, 2, 3, 4, 5, 6, 7], thus probing energy scales which are far beyond the 
reach of accelerator experiments. The simplest choice of initial conditions, namely stochastic initial 
fluctuations which are adiabatic, scalar, Gaussian, and scale-invariant, has been ruled out at the 
~3<7 level. Current observations are consistent with either a power-law initial power spectrum which 
is redder than scale invariant {n s — 1 ~ —0.04), and marginally consistent with initial conditions 
which are scale-invariant but contain contributions from tensor modes (r ~ 0.2) [8]. The next few 
years will bring a wealth of new data which will sharpen this picture considerably. 

Primordial non-Gaussianity has emerged as a particularly powerful probe of inflation due to the 
ability to rule out large qualitative classes of models. For example, there is a theorem [9, 10, 11] 
which states that in all models of single field inflation whose power spectrum is nearly scale invariant, 
the 3-point function (£(ki)£(k2)C(k3)) is observationally indistinguishable from zero in "squeezed" 
triangles (i.e. min(fcj) <C max(/cj)). However, a detectably large squeezed 3-point function is naturally 
generated in other models, such as the ekpyrotic scenario [12, 13, 14, 15]. Observational constraints 
on Gaussianity to date have mainly focused on the 3-point function and have parameterized the 
deviation from zero by three parameters f§f\ f°^f og [16, 17, 18, 19, 20, 21, 22]. The current 

WMAP constraints from [8] are: f l £f l = 32 ±21, = 26 ±140, and /^ hog = -202 ±104 (errors 

are la). 

In this paper, we will focus on local non-Gaussianity and use the notation /nl = /]vL al through- 
out. In this case, the initial curvature fluctuation is of the form £(x) = £g( x ) + |/wl(Cg( x ) 2 — 
(Cg-(x))), where Cg is a Gaussian field [23, 24, 25]. Recently, this type of non-Gaussianity has been 
studied extensively in the context of large-scale structure, beginning with a pioneering paper by 
Dalai et al. [26], which showed that scale-dependent halo bias is generated on large scales. Subse- 
quently, this prediction has been confirmed and extended, in both analytical and simulation-based 
studies [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. In particular, in [33] an improved expression 
for the Jn l dependence of the bias was introduced (the improved expression agrees with the original 
expression from [26] in the limit k — > 0). In [28], the constraint /]^£ al = 20 ± 25 was obtained from 
observations of large-scale halo clustering in SDSS (the measurement is obtained from a variety of 
tracer objects, but the statistical weight is dominated by the high-z photometric quasar sample). 
One qualitative finding which will be particularly relevant for this paper is that the scale-dependent 
bias is non-stochastic, in the sense that the correlation coefficient between halos in different mass 
bins is equal to one, after shot noise has been subtracted. 

The curvaton model is a two-field model of inflation in which the source of initial curvature 
fluctuations is not the inflaton, but a second field a whose contribution to the energy density during 
inflation is subdominant [39, 40, 41, 42, 43, 44]. Most studies of the curvaton model have only 
considered the case where the curvaton contribution to the primordial curvature fluctuation ( is 
much larger than the inflaton contribution. In this case, the curvaton field can give rise to non- 
Gaussianity of the local type at a detectable level (/nl is essentially a free parameter of the model). 
Recently, Tseliakhovich, Hirata, and Slosar considered the more general case in which the ratio £ 
of inflaton and curvaton contributions to the primordial curvature fluctuation ( can be significantly 
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different from zero [45] (see also [46] ) . Applying the same theoretical arguments which predict scale- 
dependent, non-stochastic halo bias for £ = 0, the authors argue that scale-dependent stochastic halo 
bias should be present in the more general case where Jnl and £ are both nonzero. A two-component 
hybrid model with similar observational signatures was studied in [47]; in this case £ is typically of 
order one if the non-Gaussianity is large enough to be detectable. 

The purpose of this paper is to analyze iV-body simulations whose initial conditions contain cur- 
vaton and inflaton contributions, and study the dependence of halo clustering and halo stochasticity 
on the parameters {/nl,£,} of the model. Although our main interest is the case £ / 0, we also 
present results for the curvaton model with £ = as a baseline for comparison. 

Throughout this paper we use the WMAP5+BAO+SN fiducial cosmology [48], with baryon 
density VL b h 2 = 0.0226, CDM density n c h 2 = 0.114, Hubble parameter h = 0.70, spectral index 
n s = 0.961, optical depth r = 0.080, and power-law initial curvature power spectrum k 3 P^(k)/2ir 2 = 
A^fe/fepiv)™ 3 " 1 where A 2 - = 2.42 x 10~ 9 and k piv = 0.002 Mpc^ 1 . All power spectra and transfer 
functions have been computed using CAMB [49]. 

2 Curvaton model with £ ^ 
2.1 Initial conditions 

In this subsection we review the curvaton model, in the same generality as [45] . 

The curvaton is assumed to decay before dark matter freezout, so that no dark matter isocurva- 
ture mode is generated, and the adiabatic curvature fluctuation £ is a sum of inflaton and curvaton 
contributions: 

C = Ci + Cc (i) 

We assume that the fields (i and £ c are uncorrelated and that their power spectra Pq , P^ c , are 
proportional, so that we can define a parameter £ = {Pq/Pq) 1 ^ 2 which is independent of scale. The 
power spectra of Q, £ c are thus related to the power spectrum P^ of the total curvature fluctuation 
by: 

The power spectrum P^ is taken to be of power-law form (k 3 /2TT 2 )P^(k) = /\ 2 {k/k p \ v ) ns ~ l with 
parameters A^, k p - w given in §1. 

We assume that Q is a Gaussian field, but Cc is a non-Gaussian field of "local type", i.e. 

Cc(x) = C, G (x) + p NL (l + f) 2 ( Cc 2 G (x) - <C c 2 G (x)» (4) 

where Cc,G is a Gaussian field. Non-Gaussianity of local type is generated if the curvaton potential 
V(a) is assumed quadratic in a. Throughout this paper, we will take {/tvl,c;} to be the parameters 
of the curvaton model. 
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To get some intuition for this parameterization, it is useful to note that the power spectrum, 
bispectrum, and connected higher-point functions of the initial curvature fluctuation £ depend on 
/nl and £ as follows: 

<C(ki)C(k 2 )> = Pc(fc 1 )(2vr) 3 <5 3 (k 1 +k 2 ) + 0(/2 rL ) (5) 

(C(k 1 )C(k 2 )C(k 3 )) = IfNLB^MM^^i+^ + ^ + OUltL) (6) 

5 

(C(kl)C(k 2 ) • • • C(kiv))conn = (IfNi) ^(1 + £ 2 ) "^F^ , . . . , k N )5 3 ( £ k*) + 0(/# L ) . (7) 

To lowest order in /jvx, the 3-point function is proportional to /ml, with no £ dependence 1 . This 
makes it easy to interpret the CMB bispectrum constraint from [8] as a constraint /jvl = 32 ± 21 
(la error), with no constraint on £. The CMB trispectrum constraint from [50] can similarly be 
interpreted as a constraint tnl = (1-35 ± 0.98) x 10 4 on the combination of parameters t/vl = 
(§/A r z,) 2 (l + £ 2 ), with the caveat that the bispectrum and trispectrum estimators are not statistically 
independent and there are subtleties in combining the assoicated parameter constraints [51]. 

The effect of nonzero £ is to boost the amplitude of the iV-point correlation functions (where 
N > 4) relative to the amplitude of the 3-point function. One interesting consequence is that the 
4-point function can be made large while keeping the 3-point function within observational limits 
on /tvl- (In fact, halo stochasticity, or "boosting" the amplitude of the halo-halo power spectrum 
Phh relative to the amplitude of the matter- halo power spectrum P m h, can be viewed as a formal 
consequence of boosting the primordial 4-point function relative to the 3-point function.) 

2.2 Halo clustering and the peak-background split 

The peak-background split formalism is a heuristic argument for predicting correlation functions in 
which one or more scales is large compared to the scales which are relevant for spherical collapse 
[52, 53, 54, 55, 56]. This formalism can be applied to study non-Gaussian halo clustering in the 
curvaton model [28, 45]. We will review this calculation here, and extend it by including 1-halo 
terms which will be relevant for the stochasticity results to be presented later. 

Let us write the inflaton contribution Q to the initial curvature fluctuation as a sum of long- 
wavelength and short-wavelength pieces: Q = (Qj + Q,s)- We analogously write the Gaussian field 
Cc,g as a sum (( c j + ( c ,s)- The total initial curvature fluctuation is then given by 

C(x) = C ,(x) + C c ,G(x) + ^iVL(l + £ 2 ) 2 (C c , G (x) 2 -(C c , G (x) 2 )) (8) 
= Ci,j(x) + Ccl(x) 

+Cm(x) + (l + \f NL {l + e 2 ) 2 Cc,Kx)) Ce, S (x) 

+\f N L(i + efich + c 2 s - (O + (el)) (9) 

1 We have defined /jvl in Eq. (4) with the extra factor of (l + £ 2 ) so that the 3-point function will have this property. 
The parameter f NL from [45] is related to our parameterization by f NL = /jvl(1 + £ 2 ) 2 • 
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and we have assumed the long and short wavelength parts of ( c are uncorrected. Let us make 
the approximation that the terms in the third line of (9) are negligible. (The term will not be 
important for our purposes since we will only use the small-scale component of Eq. (9); the main 
effect of the (£ a term is to change the constant part of the halo bias, which is a free parameter 
anyway.) The first line is the long- wavelength part of £, which is unchanged from the Gaussian 
case (i.e. it does not depend on /nl)- The second line is the short-wavelength part; we find that 
the effect of the non-Gaussianity is to modulate the small-scale curvaton mode £ C)S by a factor 
(1 + |/jvl(1 + £, 2 ) 2 (c,i) which depends on the long- wavelength curvaton mode £c,«- 

In the peak-background split picture, we interpret Eq. (9) as saying that the small-scale matter 
power spectrum is no longer spatially constant in a non-Gaussian cosmology, but rather a local 
quantity which varies with position. If we consider a large box at position x, then the average 
small-scale power spectrum in the box is given by 

p<w = PQ + ^+pNL(i+e) 2 CcA^) 2 P< c +0(f NL ) 

« (i + j|wi + aCc,Kx)) p c . (io) 

Following notation from [28], we will parameterize the amplitude of the small-scale power spectrum 
by ag, the RMS of the linear density field at z = with 8/i _1 Mpc tophat smoothing, and rephrase 
Eq. (10) by writing as as a function of position x: 

^8(X) = (l + pNL(l + £ 2 )Cc,i(x)) ^ • (11) 

Now let us ask how the number density of halos varies on large scales in the peak-background 
split picture. The density of halos n^(x) in a large box at position x will differ from the mean density 
fih for two reasons: first, because the local matter density p m (l + <5/(x)) in the box differs from the 
mean p m , and second because the local value of as differs from the mean via Eq. (11). Combining 
these effects we can write: 

n fc (x) = n h (l + *,(x)) (l + ¥x)^gp + pNLil + aCc,Kx)|^) (12) 

The (1 + 5i) prefactor comes from converting Lagrangian to Eulerian space. The second term is 
proportional to the (scale-independent) derivative (d log n^/ddi) of the mass function with respect 
to the background density, i.e. the Lagrangian halo bias. These two terms are present in a Gaussian 
cosmology and represent the usual halo bias which is constant on large scales. The third term repre- 
sents the effects of primordial non-Gaussianity, which gives an extra scale-dependent contribution. 
Taking the Fourier transform of Eq. (12) and dropping second-order terms, we get 

4(k) = (l + lr jW + -fc(l + 0^(c(k) 

= b G 5(k) + (1 + f)b NG (k)5 c (k) (13) 
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where 5ft (k) = rift(k)/nft is the fractional halo overdensity, and in the last line we have defined 

6G - 1 + ^ (14) 
, Vnl d log rift 

otvg(«0 = -77 — s^j ( 15 ) 

a(k, z) a log as 



a (M = on H 2 (16) 



where 

2fc 2 r(fc)£>( z) 

The quantity a(/c,z) relates the matter overdensity 5 m (k, z) to the initial curvature £(k) in linear 
perturbation theory: 5 m (k, z) = z)((k). Note that we have defined inflaton and curvaton 

contributions to the matter density by 5j(k, z) = |a(fc, z)Ci(k) and S c (k,z) = |a(/c, z)£ c (k), even 
though strictly speaking, Poisson's equation applies only to the sum of the two fields. 

The peak-background split expression (13) for 5ft (k) applies on scales which are large compared 
to scales relevant for spherical collapse. It is also incomplete, in the sense that it treats the halo 
overdensity as a continuous field, and ignores stochastic variations due to random halo locations. If 
we compute matter-halo and halo-halo power spectra using this expression, then we get: 

Pl H h(k) = [b G + b NG (k)]Piin(k) (17) 
PffiW = [(b G + b NG (k))(b' G + b' NG (k)) + ebNG(k)b' NG (k)]P iin (k) (18) 

where the primes refer to different halo masses and P\m(k) = ^a 2 (k)P^(k) is the linear theory 
matter power spectrum. We have included the superscript "2H" because these expressions omit 
1-halo terms. Using standard machinery from the halo model [57, 58, 59, 60], it is straightforward 
to calculate 1-halo contributions to these power spectra. In non-overlapping mass bins, let rii be the 
number density of halos in the i-ih bin, and let fi denote the total fraction (by mass) of dark matter 
in halos in mass bin i. Then we get: 2 

Pmm(k) = PUk) + P^m (19) 

Pmi(k) = (b < S + b% ) G (k))P ]in (k) + - (20) 

fli 



Pij(k) = (4? + b%{k))Q$ + b%{k)) + eb%{k)b%{k) PUk) + ^ (21) 

L J flj 



5, 



where we have defined 



Pnl = P^ 2 J dMM 2 n(M). (22) 

Although the 1-halo terms are generally smaller than the 2-halo terms on the angular scales we will 
study in this paper (k < 0.04 h Mpc -1 ), they are the leading souce of stochasticity (aside from 
shot noise) predicted by the halo model in the Gaussian case. In the non-Gaussian case, the term 
(^tvg^tvg) represents extra stochasticity on large scales for nonzero f^L and £ > 0. This source of 



2 These expressions neglect convolution by the halo density profiles, but this can be neglected for purposes of this 
paper, where we only study clustering on large scales (k < 0.04 h Mpc -1 ). 
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stochasticity can be understood intuitively: the non-Gaussian part of the bias traces the curvaton 
field 5 C on large scales, and if £ > 0, this field is not 100% correlated to the matter overdensity S. 

There is a useful simplification to the above expression if we assume a "universal" halo mass 
function of the form dn/dM = (p rn /M)f(v)(dv/dM). Here, v(M,z) = (5 2 c / a 2 (M , z)) , where S c = 
3(127r) 2 / 3 /20 ~ 1.69 is the threshhold for spherical collapse, and a 2 (M, z) is the variance of the 
linear matter overdensity after tophat smoothing on the scale corresponding to halo mass M. For a 
universal mass function, the relation 

dlogn <91ogn 
d log a s odi 

holds [28], so that b^G and &g are related by 

b NG (k) = -^-J N L(b G ~ 1) • (24) 

Mass functions obtained from simulations are roughly universal [61, 55, 62, 63], and we will generally 
assume that Eq. (24) holds throughout the paper. 



3 iV-body simulations 

To study halo clustering in the curvaton model, we performed collisionless TV-body simulations using 
the GADGET-2 TreePM code [64]. Simulations were done using periodic box size Rbox = 1600 

1/3 

h^ 1 Mpc, particle count N p = 1024 3 , and force softening length R s = 0.05(R^ OX /N p h ). With these 
parameters and the fiducial cosmology from §1, the particle mass is m p = 2.92 x 10 11 h^ 1 M & . 
Results in this paper were obtained from two simulations with f^L = 0, and two simulations for 
each choice of f^L G {±250, ±500} and £ G {0, 1} (for a total of 18 simulations). 

For given curvaton model parameters /atl,Cj we simulate initial conditions as follows. First, we 
simulate Gaussian fields Q and Cc,G hi Fourier space with power spectra given by Eqs. (2), (3). We 
then compute the non-Gaussian curvaton field in real space by ( c = Cc,g±!/jvl(1±£ 2 ) 2 (CcG — (Ccg))- 
(When generating the initial conditions, all Fourier transforms are computed on a grid with N p 
elements.) We apply the transfer function T{k) to the total curvature fluctuation £ = Q + ( c to 
obtain the Newtonian potential <5(fc) at the initial redshift Zi n { = 100 of the simulations. Finally, 
we obtain initial particle positions using the Zeldovich approximation [65]. (At z- m \ = 100, transient 
effects due to use of this approximation should be negligible [66].) 

We group particles into halos using an MPI parallelized implementation of the friends-of-friends 
(FOF) algorithm [67] with link length -Lfof = 0.2i?boxA"p" 1//3 . For a halo containing A^fof particles, 
we assign a halo position given by the mean of the individual particle positions, and a halo mass 
given by: 

m h = m p (A^fof - A^fof) • (25) 

The second term is recommended in [68] to minimize particle resolution artifacts when estimating 
the mass function using a FOF halo finder. 

In this paper we will analyze matter and halo power spectra, using redshifts and halo mass bins 
defined in Tab. 1. We estimate power spectra by assigning particle positions (or halo positions) to 
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a real-space grid with 1536 3 points using the Cloud-in-Cell algorithm, taking the Fourier transform, 
and averaging the power over all Fourier modes in a fc-bin. This is described in more detail in 
Appendix A. 





Mass range (h 1 M ) 


n (h 3 Mpc~ 3 ) 


/ 


b G 


z = 2 


M > 1.15 x 10 13 


3.782 x 10~ 5 


0.010 


5.222 ± 0.041 


z — 1 


1 It; v 1 nl3 ^- l\/f ^ o on i ril3 
1.10 X 1U < 1V1 < Z.oZ X 1U 

M > 2.32 x 10 13 


l.lDo X 1U 

6.346 x 10" 5 


0.038 


Z.OUZ ± U.U14 

3.470 ±0.018 


z = 0.5 


1.15 x 10 13 < M < 2.32 x 10 13 
2.32 x 10 13 < M < 4.66 x 10 13 
M > 4.66 x 10 13 


1.678 x 10~ 4 
7.585 x 10" 5 
4.498 x 10" 5 


0.035 
0.032 
0.057 


1.720 ±0.010 
2.101 ±0.014 
2.996 ±0.016 


z = 


1.15 x 10 13 < M < 2.32 x 10 13 
2.32 x 10 13 < M < 4.66 x 10 13 
4.66 x 10 13 < M < 1.02 x 10 14 
M > 1.02 x 10 14 


2.020 x 10" 4 
1.011 x 10" 4 
5.189 x 10" 5 
2.778 x 10~ 5 


0.042 
0.043 
0.045 
0.079 


1.202 ±0.008 
1.437 ±0.010 
1.783 ±0.014 
2.628 ±0.015 



Table 1: Mass bins used throughout this paper, with halo number density n, total fraction (by 
mass) / of dark matter in halos in each bin, and Gaussian bias be estimated from simulation. The 
fitting procedure used to estimate be and assign statistical errors is described in §4 and uses only 
wavenumbers k < 0.04 h Mpc -1 . 



4 Halo bias 

For a halo mass bin i, we define the bias parameter 

W*) = F^- (26) 

In this section, we will compare values of b m i(k) estimated from simulation with the predicted form: 

bmi(k) = b + ^^fNL(b - 1) • (27) 

In writing down this predicted form, we have omitted 1-halo terms derived previously in Eqs. (19)- 
(21). We find that including 1-halo terms does not qualitatively affect any conclusions from this 
section, but the simplified form in Eq. (27) is convenient for comparison with the rest of the literature. 

In Fig. 1 we show the dependence of b m i(k) on the non-Gaussianity parameters Jnl and £, for 
several choices of redshift and mass bin, on large angular scales (k < 0.04 h Mpc -1 ). 

Assigning error bars in Fig. 1 is nontrivial. On large angular scales, variations in P m i(k) and 
Pmm{k) are highly correlated (since variations in both power spectra are mainly due to sample 
variance) and therefore mostly cancel when we take the ratio to estimate b m i(k). In Appendix A we 
show in detail how to assign error bars in a way which accounts for this correlation. The resulting 
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Figure 1: Halo bias b m i(k) for selected redshifts and halo mass bins, estimated from iV-body simu- 
lations as described in Appendix A. The curves are the predicted form in Eq. (27), with fro treated 
as a free parameter which is fit from data. 
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error bars are smaller than sample variance would suggest, and result in very small statistical errors 
when fitting a given functional form of b m i(k). 

Let us separate the issue of whether the predicted form of the bias in Eq. (27) agrees with 
simulation into three separate questions. 

First, we can ask: for Gaussian initial conditions (i.e. Jnl = 0), is the bias constant on large 
scales, as predicted by Eq. (27)? If we fit for a constant bias for k < 0.04 h Mpc -1 in each redshift 
and mass bin, then we find acceptable x 2 values for each fit, indicating that the bias is indeed 
constant within the statistical errors of the simulations. This fitting procedure was used to estimate 
b G and assign statistical errors in Tab. 1. Note that the statistical error a(b G ) returned by each fit 
is typically of order ~0.01, so this \ 2 test shows that the bias is constant to percent level. 

Second, we can ask: if /jvl / but £ = 0, does the predicted form of the non-Gaussian bias 
in Eq. (27) agree with simulation? We fit for this functional form of b m i{k) over scales k < 0.04 
h Mpc -1 , treating bo as an independent free parameter for each value of /nl, i-e. we do not assume 
that the constant part of the bias bo at /nl 7^ is equal to the constant Gaussian bias be at /nl = 0. 
In [33], a slightly different fitting procedure was used: the bias b m i{k) is fit to the functional form 3 

b mi (k) = b G + -^—^f NL (b G - 1) + (A6 7 ) (28) 

where {Abj) is derived by replacing in Eq. (14) with a non-Gaussian mass function and b G is 
the Gaussian (i.e. fjy^ = 0) bias. This functional form is not precisely equivalent to Eq. (27): the 
two differ at 0(ffj L ) by the term (2S c fNL/(x(k, z))(Abi). We have not investigated whether one 
of the functional forms is a better fit than the other, but we anticipate that an 0(ffj L ) difference 
will be negligible for values of /nl which are observationally relevant. However, we do find that 
including an O(Jnl) term which is constant in k, either via the last term in Eq. (28) or by allowing 
6o to differ from b G in Eq. (27), is needed to obtain a good fit to simulation. The magnitude of this 
scale-independent Jn l correction we find is in qualitative agreement with what one would get using 
the non-Gaussian mass-function of [69] in Eq. (14). 

We find that for some choices of redshift and halo mass bin, the fits return bad x 2 values (Tab. 2), 
i.e. we find statistically significant disagreement between the simulations and the predicted form of 
the bias in Eq. (27). Detailed inspection of the bad fits shows that, in all cases, the prediction 
tends to overestimate the magnitude of the non-Gaussian bias on large scales, but the discrepancy 
between simulation and prediction is only ~10% of the total non-Gaussian bias in the worst case. 
The theoretical assumptions made in deriving the non-Gaussian bias, i.e. a universal mass function 
and the peak-background split relation between halo bias and the mass function, also fail at roughly 
this level [70, 71], so a «10% discrepancy is not at all surprising. 

Third, we can ask whether the bias b m i(k) is independent of £ for fixed /atl, as predicted by 
Eq. (27). To test this, we let b(k) and b'(k) denote estimates of the bias from two ./V-body simulations 

3 Note that the expression for Ab in Eq. (9) of [33] also includes a term denoted (b(M)/3 m (k, /nl)) which corresponds 
to /jv_l dependence of the matter power spectrum P mm (fe). We do not include this term because we define the bias 
to be b m i(k, /nl) = P m i{k, fNL)/Pmm{k, /nl) with an /jvL-dependent denominator. Such a term would be needed if 
the bias were defined as (P m i(k, f NL )/ P mm (k,0)) or as (P m i(k, f NL )/Piin(k)). 
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f NL = 500 


Inl — — 


500 




Mass range (/i 1 M ) 


bo 




bo 




z = 2 


M > 1.15 x 10 13 


4.228 ± 0.020 


61.6/13 


6.560 ±0.083 


70.7/13 


z = 1 


1.15 x 10 13 < M < 2.32 x 10 13 


2.238 ±0.008 


29.6/13 


2.878 ±0.023 


65.9/13 




M > 2.32 x 10 13 


3.012 ±0.011 


42.9/13 


4.141 ±0.033 


46.4/13 


z = 0.5 


1.15 x 10 13 < M < 2.32 x 10 13 


1.623 ±0.007 


28.3/13 


1.867 ±0.015 


40.4/13 




2.32 x 10 13 < M < 4.66 x 10 13 


1.940 ±0.010 


11.8/13 


2.346 ±0.022 


25.4/13 




M > 4.66 x 10 13 


2.666 ±0.010 


36.4/13 


3.448 ± 0.026 


31.0/13 


z = 


1.15 x 10 13 < M < 2.32 x 10 13 


1.193 ±0.006 


27.6/13 


1.229 ±0.010 


17.3/13 




2.32 x 10 13 < M < 4.66 x 10 13 


1.386 ± 0.008 


7.7/13 


1.524 ±0.014 


23.2/13 




4.66 x 10 13 < M < 1.02 x 10 14 


1.675 ±0.010 


14.9/13 


1.907 ±0.018 


30.5/13 




M > 1.02 x 10 14 


2.403 ±0.010 


23.8/13 


2.946 ±0.022 


20.1/13 



Table 2: Best-fit values of bo, and \ 2 values for the fit, when fitting the predicted form of the non- 
Gaussian bias in Eq. (27) to estimates of the bias b m i(k) from simulation, with error bars assigned 
as described in Appendix A. A few of the fits return bad \ 2 values; in these cases we find that 
Eq. (27) overpredicts the non-Gaussian bias by «10%. 

with £ = and £ = 1 (and the same value of Jnl)- We define a x 2 statistic by summing over fc-bins: 

y2 = V m - h ' {k)? (29) 

X ^ Var(A6(fc))±Var(A6'(A:)) 1 ; 

For almost all redshifts and mass bins, we find acceptable \ 2 values, indicating that the bias estimates 
from the two simulations are consistent within statistical errors. (There is one exception: we find 
an anomalous \ 2 f° r Inl = —500 and z = 2, but inspection of the bad fit shows that the bias only 
differs by ps 10% between £ = and £ = 1. The bias is larger in the £ = 1 case, at low k.) 

Our conclusion in this section is that the prediction for the non-Gaussian halo bias in Eq. (27) is 
an impressive fit to the simulations across a wide range of curvaton model parameters. Although we 
do detect statistically significant deviations from the prediction at the « 10% level, this is typical 
for results based on general arguments such as the peak-background split. 



5 Halo stochasticity 

For halo mass bins i, j, we define the stochasticity parameter 

s _ Pjjjk) — Sjj/rii Pmi{k)P m j{k) 

ij{ ' ~ p Ik) p jw ■ [ ' 

If the halos are perfectly non-stochastic tracers of the dark matter (or more precisely, if the only 
source of stochasticity is shot noise) then both the diagonal (i.e. i = j) and non-diagonal (i.e. i ^ j) 
components of rij(k) will be zero. As discussed in §2.2, the halo model predicts that the leading 
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contribution to rtj(k) on large scales arises from 1-halo terms in the matter-halo power spectra. If 
£ > 0, then we expect large since the non-Gaussian part of the bias will be stochastic. 

In Fig. 2, we show estimates of the diagonal components ra(k) from simulation, for several 
choices of redshift, mass bin, and model parameters (/jvx> £)• Similarly to the case of halo bias from 
the previous section, assigning error bars is nontrivial, since estimates of P mm , Pmi, and Pij are all 
highly correlated on large scales. Our procedure for estimating the stochasticity and assigning error 
bars is given in Appendix A and results in error bars which are much smaller than sample variance 
would suggest. 

In this section, we will compare the stochasticity estimated from simulations to the predicted 
form 



which follows from the halo model calculations in §2.2. As in the previous section, we will separate 
the issue of whether this prediction agrees with simulation into three separate questions. 

First, does the prediction in Eq. (31) agree with simulations in the Gaussian case (/nl = 0)? 
Surprisingly, we do not even find agreement at a qualitative level. Although there are a few choices 
of redshift and mass bin where the halo model prediction roughly fits the data, there are more cases 
where there is no resemblance (example good and bad fits are shown in Fig. 3). 4 

In this paper, we have not attempted to propose a general model for r„ in the Gaussian case. 
Therefore, when we compare the prediction for rij to simulation in the non-Gaussian case, our 
approach is to estimate the change in stochasticity Arjj between non-Gaussian and Gaussian initial 
conditions, as a function of (/nl, £)> an d compare with the prediction for Arij obtained from Eq. (31). 

The second question we can ask is, how does the stochasticity ru depend on /jvl, in the case 
£ = 0? In this case, we find that Eq. (31) predicits an /jvl dependence which is small compared to 
the statistical errors of our simulations. To test this prediction, we define a \ 2 statistic by: 



where fu(k) and r'u{k) denote stochasticity estimates from A-body simulations with /nl = and 
/nl 7^ respectively (taking £ = in the non-Gaussian simulation). For /nl = ±500 and all 
redshifts and mass bins considered in this paper, we find good \ 2 values, i.e. no statistically significant 
dependence of the stochasticity ru on /jvl (provided £ = 0). This agrees with the halo model 

4 A recent paper [72] also compared stochasticity predictions in the halo model with iV-body simulations in the 
Gaussian case. It was found that the smallest eigenvalue of the stochasticity matrix nj is predicted accurately by 
the halo model. (The smallest eigenvalue is particularly relevant since it corresponds to a halo mass weighting with 
reduced shot noise [73].) However, it can also be seen (Fig. 12 of [72]) that the halo model does not accuractely predict 
the largest eigenvalue, which implies that some matrix elements m are not accurately predicted. In this paper, we 
have concentrated on the diagonal ru and do not generally find good agreement with the halo model. 




(31) 




(32) 
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Figure 2: Diagonal components ra of the stochasticity statistic defined in Eq. (30), estimated from 
A-body simulations as described in Appendix A, for varying choices of redshift, halo mass range, 
and curvaton model parameters (Jnl, £)• 
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Figure 3: Stochasticity parameter ra estimated from Gaussian simulations (error bars), with halo 
model prediction shown for comparison (curves). In general, we do not find that the halo model 
accurately predicts ra. An example of a redshift and mass bin where the halo prediction disagrees 
with simulation (z = 2 and M > 1.15 x 10 13 }i~ 1 Mq) and an example where the two agree [z = 0.5 
and M > 4.66 x 10 13 h^Mp,) are shown. 
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prediction (31), even though the halo model does not correctly predict the actual value of ru as 
previously remarked. 

Third, we can ask, how does the stochasticity depend on £? Since we do not have a model for 
the Gaussian stochasticity, we ask whether the quantity 



i.e. the excess stochasticity over Gaussian, is correctly modeled by the peak-background split pre- 
diction: 



(This is actually an approximation to the prediction obtained by differencing Eq. (31) between values 
of (/jvLjOj Du t we find that this approximation is within statistical errors of the simulation, so it is 
a convenient simplification.) 

We find that the prediction in Eq. (34) systematically overestimates (Aru)- Empirically, we find 
that if we scale the prediction for (Aru) by a multiplicative constant q, then the modified prediction 



is an excellent fit to the simulations (i.e. when q is treated as a free parameter, all fits have good 
X 2 values), for all choices of redshfit, halo mass bin, and /jvl- In Fig- 4, we show an example fit; it 
is seen that the peak-background split overpredicts the level of stochasiticity, but an excellent fit is 
obtained by simply scaling the peak-background split prediction. In Tab. 3, we tabulate values of q 
obtained by fitting Eq. (35) to the simulations, together with statistical errors from the fits. 

One general trend evident in this table is that q is an increasing function of /nl- This simply 
means that we are considering large enough values of Jnl that we are sensitive to 0(fff L ) terms, if 
we think of (Aru) as a power series in /jvl- However, it is clear from Tab. 3 that if we interpolate to 
Inl — > 0, most values of q are still < 1. This indicates that, even to leading order 0(f^ L ), we are 
seeing ~ 30% disagreement between the prediction (34) and simulation. We have not attempted to 
study 0(ffj L ) terms in detail, or compare them between theory and simulation, since we are already 
seeing disagreement at leading order 0(f^ L ). 

Since the typical value is q « 0.7, our interpretation is that the peak-background split generally 
overpredicts non-Gaussian stochasticity by ~30%. Although some discrepancy is expected (e.g. in 
the previous section we found < 10% discrepancies in predictions for halo bias), this level of disagree- 
ment is somewhat uncomfortable and should probably be incorporated when constraining £ from 
observations. It is not clear how to interpret this discrepancy theoretically, or whether it is related 
to the discrepancy that we found previously in the Gaussian case (Fig. 3). Modeling stochasticity 
using a semianalytic framework such as the peak-background split or halo model appears to be more 
difficult than modeling halo bias. 

6 Discussion 

In this paper, we have compared semianalytic predictions for halo clustering to iV-body simulations, 
in the two-field inflationary model from [45], in which the initial curvature fluctuation is a sum of 



Ar. 



■a = ru(k, /nl,0 - ru(k, f NL = 0) 



(33) 




(34) 




(35) 
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Figure 4: Change in stochasticity parameter Aru = ra(k, Snl,0 — r u(k, Snl = 0) between the 
curvaton model with (Snl,0 = (500,1) and the Gaussian case, estimated from iV-body simula- 
tions. The peak-background split (solid curve) overpredicts the level of stochasticity, but excellent 
agreement with simulation is obtained by scaling the prediction by q = 0.42 (dotted). When the 
parameters (Snl,£,z) and the halo mass range are varied, we find that scaling the peak-background 
split prediction always provides a good fit, but the value of q varies, as shown in Tab. 3. 





Mass range (h 1 M ) 


Snl = 500 


Snl = 250 


Snl = -250 


Snl = -500 


z = 2 


M > 1.15 x 10 13 


0.98 ±0.07 


0.88 ±0.08 


0.62 ±0.06 


0.42 ±0.03 


z = 1 


1.15 x 10 13 < M < 2.32 x 10 13 
M > 2.32 x 10 13 


0.79 ±0.09 
0.83 ±0.07 


0.83 ±0.12 
0.70 ±0.08 


0.67 ±0.09 
0.66 ±0.07 


0.46 ± 0.04 
0.51 ±0.04 


z = 0.5 


1.15 x 10 13 < M < 2.32 x 10 13 
2.32 x 10 13 < M < 4.66 x 10 13 
M > 4.66 x 10 13 


1.01 ±0.18 
0.80 ±0.15 
0.81 ±0.09 


0.92 ±0.29 
0.58 ±0.22 
0.79 ±0.12 


0.45 ±0.19 
0.73 ±0.19 
0.80 ±0.10 


0.57 ±0.10 
0.48 ±0.08 
0.51 ±0.05 


z = 


1.15 x 10 13 < M < 2.32 x 10 13 
2.32 x 10 13 < M < 4.66 x 10 13 
4.66 x 10 13 < M < 1.02 x 10 14 
M > 1.02 x 10 14 


1.37 ±0.80 
1.35 ±0.44 
0.71 ±0.26 
0.79 ±0.13 


1.06 ± 1.12 
1.57 ±0.77 
0.90 ±0.49 
0.93 ±0.21 


1.00 ±1.41 
0.82 ±0.59 
1.12 ±0.41 
0.73 ±0.15 


0.90 ±0.51 
0.58 ±0.25 
0.63 ±0.17 
0.53 ±0.07 



Table 3: Values of the g-parameter, defined in Eq. (35), obtained from iV-body simulations for 
various values of Snl, redshift, and mass bin. (We take £ = 1 throughout) 



Gaussian and non-Gaussian contributions. This model is parameterized by Snl, which corresponds 
to the amplitude of the 3-point function in squeezed triangles, and a parameter £ which corresponds 
to the ratio of inflaton to curvaton fluctuations and boosts the 4-point and higher functions relative 



15 



to the 3-point function. Note that curvaton models also generally predict a cubic contribution of 
the form (qnlCg) t° the initial curvature. We have not considered such a term here since the peak- 
background split analysis differs significantly from the Jnl case, and defer study of the g^L term to 
future work [74]. 

The halo bias b(k) = P m h(k) / P mm (k) in simulation is found to agree very well with the peak- 
background split prediction (27) on scales k < 0.04 h Mpc -1 , for a range of redshifts and halo 
masses. In the Gaussian case, the bias is constant in k at the percent level. In the non-Gaussian 
case, we find deviations from the functional form for b(k) predicted by the peak-background split 
which are small (< 10%) but statistically significant for our simulation volume. We interpret this 
as agreement with the prediction, since the peak-background split is expected to break down at the 
~10% level. 

We also compare the shot noise subtracted halo stochasticity parameter 



measured in simulation to semianalytic predictions. In this case the results are more puzzling; some 
of the semianalytic predictions are confirmed and others are not. In the Gaussian case (Jnl = 0), 
the halo model makes a prediction for r(k) (the leading contribution is from the 1-halo term in P m h), 
but this prediction does not consistently agree with simulation. The excess stochasticity (relative to 
the halo model prediction (31)) observed in simulations is consistent with an additive Poisson-likc 
contribution to the halo-halo power spectrum (i.e. APhh(k) is independent of k on large scales, but 
can depend on redshift and halo mass). 

In the non-Gaussian case, the halo model also predicts that the stochasticity does not depend 
on Inl (if £ = 0); we find that this is true in the simulations. The last prediction is a specific 
functional form (34) for the difference stochasticity (Ar„) in stochasticity between a model with 
£ > and a Gaussian cosmology. We find that the simulations deviate from this prediction by an 
overall multiplicative factor q ~ 0.7, a significant enough disagreement that it should probably be 
incorporated when constraining the two-parameter curvaton model from observations. It is not clear 
whether the disagreements between theory and simulation in the Gaussian and non-Gaussian cases 
are related; it is not straightforward to compare the two since the excess stochasticity appears to be 
"additive" in the Gaussian case and "multiplicative" in the non-Gaussian case. 

Large-scale halo clustering has emerged in the last few years as one of the most powerful probes of 
primordial non-Gaussianity of one of the local types (i.e. either /jy^ 11 , S^zf 1 ' or t ne two-field local type 
considered here). In this paper, we have confirmed qualitative predictions from the peak-background 
split ansatz, showing that the peak-background split picture is very useful for relating local- type 
primordial non-Gaussianity to observations. However, detailed comparison reveals differences which 
are large enough (« 30% in this case) to be important for data analysis, highlighting the need 
for simulations. The next few years should bring a mixture of theoretical, simulation-based, and 
observational work which will greatly sharpen our observational constraints on the physics of the 
early universe. 




2 
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A Estimators for b and r 

Throughout this paper, we have given estimates of the halo bias bi(k) = P m i(k) /P mm (k) and stochas- 
ticity ru(k) = (Pu(k) — I / rii) / P mm (k) — (P m i(k) / P mm (k)) 2 , with error bars that are used for param- 
eter fitting and computing x 2 statistics. In this appendix we describe our estimator methodology, 
in particular the calculation of error bars. 

In principle, error bars could be assigned by running multiple iV-body simulations and using the 
Monte Carlo scatter in estimates of b or r, but this is impractical since the "error on the error" 
would be y/2/N mc where N mc is the number of iV-body simulations, so computing error bars with 
10% accuracy would require running N mc 200 simulations. Therefore, an analytic prescription for 
the error bars is necessary. 

One can see intuitively that in the limit of zero stochasiticity (r — > 0) and zero shot noise 
(rii —7- oo), the statistical errors on b and r should go to zero, since there will be no scatter around 
the mean relation Pa = biP m i = bfP m m- Put another way, the statistical errors should not receive 
contributions from sample variance, since sample variance cancels when we take ratios of power 
spectra. This will be reflected in our final expressions for the error bars (Eqs. (50) and (51) below), 
in which all terms contain prefactors of r or (1/n). 

In a finite volume V we use the Fourier conventions 

<5(k) = J d 3 x<5(x)e ik ' x (37) 
<5(x) = F" 1 ^(kje-** (38) 

k 

With these conventions, the infinite-volume two-point function (5(k)<5(k')*) = P(/c)(27r) 3 (5 3 (k — k') 
becomes 

(S(k)S(k')*) = VP{k)5w (39) 
In a fc-bin b, we estimate power spectra using the estimator 

k ket 

where = X^kefe ^ ^ s ^ ne numrj er of Fourier modes in the fc-bin, and indices a, (3 can denote either 
the matter overdensity field or a halo mass bin i. 

For the rest of the appendix, we fix the redshift, the fc-bin b, and the halo mass bin i. We will 
make the approximation that variation in power spectra across a single fc-bin is not important, and 
use the compressed notation P a/ 3 = P a p(k). We will also use compressed notations n = rii for the 
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halo number density, b = P m i/P mm for the bias, and r = (Pa - l/n,)/P mm - {P m i/P mm ) 2 for the 
stochasticity. 

Define estimators for halo bias and stochasticity by: 



p . 

1 mi 
p 

1 mm 



r = 



N k~ 2 \ p. _ {N k -2\ 1 



N k -\J \ N k J n 



N k 2 / Pmi 



*k-l \P„ 



(41) 
(42) 



The factors of (N k — 2)/(N k — 1) and (N k — 2)/N k in the second equation are ad hoc for now, but we 
will show (Eqs. (48), (49) below) that they ensure that f is an unbiased estimator of r in the case 
where N k is not » 1. 

To calculate Var(fo) and Var(f), we will make the approximation that the matter overdensity 5 m 
and halo overdensity Si are Gaussian fields. 5 We would first like to characterize the joint PDF of 
the power spectrum estimators P mm , Pmi and Pa. Let us first define normalized fields 

^1 = Pmm $m 

( iV 1/2 

62 = \rP mm + -\ (Si-bS m ) (43) 

whose power spectra are normalized to P\\ = P\2 = 1 and P\2 = 0. The two sets of power spectrum 
estimators are related by: 



Pmm — PmmP\\ 



12 ( 1 V /2 " 

Pmi = bP mm P\\ ~\~ Pmm [ r Pmm ~\~ ~ J P\2 

Pa = b 2 P mm Pn + 2bPUl [rP mm + ' A2 + (rP mm + P22 



(44) 



The joint PDF of the random variables (N k Pn), (-/V fe P 12 ), and (N k P22) is well-studied; it is known 
as the Wishart distribution. A convenient way to characterize this distribution is via the Bartlett- 
Cholesky decomposition, which states that: 



N k Pi2 N k P 22 




(45) 



where the random variable xo is x 2_ distributed with Nk degrees of freedom, xi 1S X 2 -distributed 
with (N k — 1) degrees of freedom, a is distributed as a unit Gaussian, and the three variables xoi 
Xi, and a are statistically independent. 



5 This assumption may appear inconsistent, since are considering non-Gaussian initial conditions, and in addition 5, 
will be a non-Gaussian field even if the initial conditions are Gaussian. However, we will only use the final expressions 
for Var(6) and Var(f) on large scales (k < 0.04 h Mpc -1 ), and on these scales the fields S m , Si should be approximately 
Gaussian. 
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Combining Eqs. (44) and (45), we find that 



Pr, 



PmmXO 








N k 








PmmXO 


b + 
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N k 


n Pmm 


PmmXO 


b + 


+ 


1 


N k 


n Pmm 



1/2 



ax 



1/2 



ax 



1/2 



1/2 



. PmmXl ( ' 
H 77- — r + 



iV fc 



(46) 



These equations, combined with the statement at the end of the previous paragraph which gives the 
joint PDF of xoj Xi> an d a, completely characterize the sampling PDF of the estimators P mm , P m i 
and Pa. 

Armed with this characterization, it is easy to calculate Var(6) and Var(f). We write the esti- 
mators b and r (defined in Eqs. (41), (42)) in terms of the variables xo, Xi, and a (using Eq. (46)), 
obtaining: 



b = b+ [r + 



nP„ 



1/2 



r = 



JV fc -l 



r + 



nP„ 



ax 

Xi 
Xo 



1/2 



N k 



nP„ 



1 

Xo 



(47) 



It is then straightforward to calculate the mean and variance of the quantities on the right-hand 
side, obtaining: 



(b) 
(f) 

Var(6) 
Var(f) 



b 

r 



N k -2 
2 

N k -A 



r + 



{r 2 ) + 



(48) 
(49) 

(50) 

(N k - l)(N k - 4) V r + nP mm ) (51) 
2 m r(m+N/2)/T(N/2), where x is x 2 -distributed 



nP, 



mm , 

2(N k 



For these calculations the expectation value (x m ) 
with N degrees of freedom, is useful. 

This calculation is exact even in the case where Nk is not 3> 1 , and this level of precision appears 
to be necessary, e.g. we find that if terms of order (l/N k ) in the variance are neglected, then a few 
jackknife tests in Tab. 4 below fail. We also note that both estimators are unbiased (i.e. (b) = b and 
(f) = r, justifying the factors of (Nk — 1) and (Nk — 2) in the definition (42) of f, which were ad hoc 
until now. (The estimators would be biased if these factors were omitted.) 

There is one final wrinkle: in order to apply these expressions for Var(6) and Var(f), we need to 
know the stochasticity r. In this paper we do not propose a general model for stochasticity, finding 
for example that the nonzero stochasticity seen in simulations for Gaussian initial conditions is not 
fit well by the halo model prediction (§5), and makes a non-negligible contribution to the estimator 
variance. Therefore, we infer the stochasticity directly from the simulation itself, by making the 
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approximation 



Pa - 1/n 



(52) 



on the right-hand sides of Eqs. (50), (51). 

Putting the results of this appendix together, our final estimates for Var(6) and Var(f) are given 

by: 



Var(6) 
Var(f) 



1 



N k -2 
2 

N k -A 



P 



p 



P / VP 

1 mm I X 1 mm 



(53) 



1/n 



+ 



2(N k - 2) 
(N k -l)(N k -A) 



Pa 



We have used these expressions to assign error bars throughout this paper, e.g. in Figs. 1 and 2. 

We conclude this appendix with an end-to-end test of our estimates for Var(fe) and Var(f). If we 
run two independent iV-body simulations with the same values of fpjL and £, then the differences 
(b — b ) between bias estimates, and differences (f — r') between stochasticity estimates, should be 
consistent with zero. This jackknife test is intended to check correctness of the error bars without 
requiring a model for the expected values of the bias and stochasticity. 

More precisely, for each halo mass bin i, we define a x 2 statistic by summing over fc-bins b: 



X 



ru(b)-r'u(b) 



Var(f ii (6)) + Vax(f; i (6)) 



(54) 



and analogously with the stochasticity estimator r replaced by the bias estimator b. Results from 
the jackknife tests are shown in Tab. 4. The results appear consistent with % 2 statistics, indicating 
that our error estimates are accurate. (The most anomalous x 2 value in the table is 27.6 with 14 
degrees of freedom, corresponding to a p-value of 1.6%. Since there are 100 entries in the table, an 
entry which is anomalous at this level is expected.) 
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X 2 for b m i 


jackknife test (AW = 


14) 






Mass range (/i 1 M ) 


Inl = 


fNL = 500 


!nl = -500 


f NL = 500 


fNL = -500 






C = o 


£ = o 


C = o 


e = i 




£ = 1 


z = 2 


M > 1.15 x 10 13 


23.4 


27.6 


16.7 


14.3 


11.1 


z = 1 


1.15 x 10 13 < M < 2.32 x 10 13 


16.5 


13.9 


7.4 


16.1 




15.1 




M > 2.32 x 10 13 


16.0 


14.4 


24.1 


13.1 




6.9 


z = 0.5 


1.15 x 10 13 < M < 2.32 x 10 13 


17.1 


19.1 


22.2 


10.1 




9.8 




2.32 x 10 13 < M < 4.66 x 10 13 


9.4 


7.0 


9.4 


19.5 




20.9 




M > 4.66 x 10 13 


19.3 


24.0 


12.3 


8.7 




7.9 


z = 


1.15 x 10 13 < M < 2.32 x 10 13 


16.7 


13.1 


11.4 


14.2 




22.5 




2.32 x 10 13 < M < 4.66 x 10 13 


25.2 


14.5 


16.5 


20.0 




8.7 




4.66 x 10 13 < M < 1.02 x 10 14 


12.3 


14.5 


11.4 


13.6 




12.1 




M > 1.02 x 10 14 


21.6 


10.9 


18.1 


9.5 




12.1 









X 2 for r a 


jackknife test (iVdof = 14) 






Mass range (h l M & ) 


fNL = 


fNL = 500 


f NL = -500 


fNL = 500 


f N L = -500 






e = o 


£ = o 


e = o 


e = i 




£ = i 


z = 2 


M > 1.15 x 10 13 


16.9 


14.0 


19.5 


9.6 


7.9 


z = 1 


1.15 x 10 13 < M < 2.32 x 10 13 


12.7 


10.7 


6.8 


10.9 




10.7 




M > 2.32 x 10 13 


7.3 


10.2 


13.2 


9.3 




8.8 


z = 0.5 


1.15 x 10 13 < M < 2.32 x 10 13 


16.8 


11.8 


9.0 


8.2 




18.9 




2.32 x 10 13 < M < 4.66 x 10 13 


9.4 


9.3 


5.2 


14.2 




21.8 




M > 4.66 x 10 13 


6.7 


12.1 


8.3 


12.0 




9.7 


z = 


1.15 x 10 13 < M < 2.32 x 10 13 


9.8 


10.8 


15.7 


8.3 




9.9 




2.32 x 10 13 < M < 4.66 x 10 13 


7.2 


13.3 


9.0 


12.0 




7.1 




4.66 x 10 13 < M < 1.02 x 10 14 


14.9 


7.8 


9.8 


9.7 




7.6 




M > 1.02 x 10 14 


24.4 


20.0 


18.3 


15.9 




21.9 



Table 4: Jackknife tests for the bias estimator b m i and stochasticity estimator fjj. The x 2 values are 
consistent with statistical expectations, showing that the variances Var(6 m j) and Var(fjj) have been 
correctly estimated. 
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